rm(list=ls())

# Parameters
Overline.v <- 3
Mu <- 2
Sd <- 1
Cost <- 1.5
W.1 <- 1
Boost <- 0.1
Decline <- 1
Underline.v <- ((W.1 + Boost) - (W.1 - Cost))/Decline

Shape1.lib <- 2
Shape2.lib <- 5

# Functions
right.hand.side.eq <- function(v, 
                               cost = Cost, w.1 = W.1, 
                               boost = Boost, decline = Decline){
  w.0 <- w.1 + boost - decline*v
  lhs <- cost / (w.1 - w.0)
  return(lhs)
}

left.hand.side.eq <- function(hat.v, 
                              shape1, shape2, 
                              overline.v = Overline.v,
                              mu = Mu, sd = Sd){
  posterior <- (1-pnorm(overline.v, mean = mu, sd = sd)) / 
    (1-pnorm(hat.v, mean = mu, sd = sd))
  out <- pbeta(q = posterior, 
                       shape1 = shape1, shape2 = shape2)
  return(out)
}

equ.condition <- function(v, shape1, shape2,
                          cost = Cost, w.1 = W.1, 
                          boost = Boost, decline = Decline, 
                          underline.v = Underline.v, 
                          overline.v = Overline.v,
                          mu = Mu, sd = Sd){
    rhs <- right.hand.side.eq(v, 
                             cost = cost, w.1 = w.1, 
                             boost = boost, decline = decline) 
    lhs <- left.hand.side.eq(v,
                             shape1 = shape1, shape2 = shape2, 
                             overline.v = overline.v,
                             mu = mu, sd = sd)
    out <- lhs - rhs
}

equ.threshold <- function(shape1, shape2, cost = Cost, w.1 = W.1, 
                          boost = Boost, decline = Decline, 
                          underline.v = Underline.v, 
                          overline.v = Overline.v,
                          mu = Mu, sd = Sd){
  out <- uniroot(equ.condition, 
                 lower = underline.v, 
                 upper = overline.v, 
                 shape1 = shape1, shape2 = shape2,
                 cost = cost, w.1 = w.1, 
                 boost = boost, decline = decline, 
                 underline.v = underline.v, 
                 overline.v = overline.v,
                 mu = mu, sd = sd)$root
  return(out)
}

# Figure 
v.full <- seq(1, 4.3, length.out = 1000)
v.restricted <- seq(Underline.v, Overline.v, length.out = 1000)
labels.x.type <- c(expression(underline(v)), expression(paste(widehat(v), '*')), expression(bar(v)))

pdf('./figures/fig_1b.pdf', width = 7, height = 7)
plot(v.restricted, left.hand.side.eq(hat.v = v.restricted, 
                                     shape1 = Shape1.lib, 
                                     shape2 = Shape2.lib), 
     type = 'l', lwd = 4, ylim = c(0, 1), col = 'grey',
     xaxt = 'n',
     main = '', 
     xlab = 'Threshold for applying', ylab = 'Probability', 
     cex.main = 1.5, cex.lab = 1.5)
axis(side = 1, at = equ.threshold(shape1 = Shape1.lib, shape2 = Shape2.lib),
     labels = expression(paste(widehat(v), '*')), 
     cex = 1.5)
lines(v.restricted, right.hand.side.eq(v = v.restricted),
      lwd = 4, col = 'black')
segments(x0 = equ.threshold(shape1 = Shape1.lib, shape2 = Shape2.lib), 
         y0 = 0, 
         x1 = equ.threshold(shape1 = Shape1.lib, shape2 = Shape2.lib), 
         y1 = right.hand.side.eq(v = equ.threshold(shape1 = Shape1.lib, shape2 = Shape2.lib)), 
         lty = 2, lwd = 3)
legend(x = 2.18, y = 0.25, 
       col = c('black', 'grey'), 
       lwd = c(3, 3), bty = 'n',
       legend = c(expression(paste('Stakes-adjusted costs: ', frac(c,w[1]-w[0](hat(v))))),
                  expression(paste('Probability accept: ', 
                                   H(q(hat(v)))))), 
       cex = 1)
dev.off()